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Abstract 

A contemporary and exciting application of Grobner bases is their use in computational biology, particularly in 
the reverse engineering of gene regulatory networks from experimental data. In this setting, the data are typically 
limited to tens of points, while the number of genes or variables is potentially in the thousands. As such data 
sets vastly underdetermine the biological network, many models may fit the same data and reverse engineering 
programs often require the use of methods for choosing parsimonious models. Grobner bases have recently been 
employed as a selection tool for polynomial dynamical systems that are characterized by maps in a vector space 
over a finite field. 

While there are numerous existing algorithms to compute Grobner bases, to date none has been specifically 
designed to cope with large numbers of variables and few distinct data points. In this paper, we present an 
algorithm for computing Grobner bases of zero-dimensional ideals that is optimized for the case when the number m 
of points is much smaller than the number n of indeterminates. The algorithm identifies those variables that are 
essential, that is, in the support of the standard monomials associated to a polynomial ideal, and computes the 
relations in the Grobner basis in terms of these variables. When n is much larger than m, the complexity is 
dominated by nm 3 . The algorithm has been implemented and tested in the computer algebra system Macaulay 2. 
We provide a comparison of its performance to the Buchberger-Moller algorithm, as built into the system. 

Keywords: Grobner bases, Buchberger-Moller algorithm, essential variables, run-time complexity, computational 
biology applications 

1 Introduction 

The theory of Grobner bases has been an active field of study in the last four decades, beginning with the seminal 
work of Buchberger [6] . A problem of particular interest has been the development of algorithms for computing 
Grobner bases. The first algorithm, proposed by Buchberger, has time complexity that is doubly exponential in 
the number of variables [TJ. Since then, several improvements to Buchberger's algorithm have been proposed, as 
well as a number of alternative methods for certain classes of ideals. 

Many of the improvements focus on two aspects. The first is coefficient growth when computing Grobner bases 
in a field of characteristic (for example, see [5]). The second is Buchberger's Criterion, which states that 

Q 

"A set G = {<7i, . . . , g r } C I is a Grobner basis for I if and only if the S-polynomial S(gi,gj) is for 
all 1 < i,j < r." 

The Optimized Buchberger Algorithm [9] proposed by Caboara et al. and Faugere's F4 and F5 [131 114] are 
instances of methods that seek to minimize the number of S'-polynomials to be computed. While they still have 
exponential complexity in the worst case, in practice their performance renders them efficient alternatives to the 
original Buchberger algorithm. 

For zero-dimensional ideals, several methods have been described and implemented. In [8], the authors pre- 
sented the Buchberger-Moller algorithm (BMA) for computing the reduced Grobner basis for the vanishing ideal 
of a variety V over a field. This algorithm eliminates the need to compute S'-polynomials and instead performs 
Gaussian elimination on a generalized Vandermonde matrix. Its complexity is quadratic in the number of vari- 
ables and cubic in the number of points in V (for details, see [20, 21, 22]). It has been implemented in publicly 
available computer algebra systems such as CbCcA [10] and Macaulay 2 [15|- The BMA was later generalized to 
noncommutative rings [3] ■ Abbott et al. [I] described a modular version of the BMA for the case when fc = Q. 
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There are other algorithms for zero-dimensional ideals which have been developed for particular settings. Farr 
and Gao presented an algorithm based on a generalization of Newton interpolation in |12] . While the complexity 
is exponential in the number n of variables, the algorithm has been optimized for the case in which n is small as 
compared to the number of points. Lederer proposed a method for lexicographic term orders which gives insight 
into the structure of the Grobner basis [19] . 

A recent and exciting development in the theory of Grobner bases is their use in computational biology. 
For instance they have been used in the identification of critical points of maximum likelihood functions in 
phylogenetic-tree reconstruction |16| . Grobner bases have also been employed as a selection tool for polyno- 
mial dynamical systems (PDSs) in the study of gene regulatory networks |18| and protein signal transduction 
networks [3]. 

In applications to molecular biology, networks often consist of n biochemicals, such as gene products or metabo- 
lites, with changing concentration levels. In [18] a method was proposed to reverse engineer biochemical networks, 
where the levels are mapped to a finite field k — ¥ p for some prime p. In this setting, networks are modeled as 
PDSs, which generalize the widely studied Boolean networks (see [T7] for an introduction). Concentration levels 
are recorded in a vector in k n , and the data consists of input-output pairs (si,Si+i) £ k n x k n , where Sj is a 
vector describing the state of the network at time i, for i = 1, . . . , m. The input vectors can be viewed as an affine 
variety V C k n , and a family of models represented as PDSs is constructed in terms of the vanishing ideal of V. 
Grobner bases are then used to select the most parsimonious PDS from this collection. In these applications, the 
number n is typically in the hundreds to thousands, whereas the number m is at best on the order of tens of 
measurements. 

Below we describe an algorithm for computing Grobner bases for zero-dimensional ideals (i.e., vanishing ideals) 
in a polynomial ring R. This algorithm is specialized for the case when the number m of distinct points is much 
smaller than the number n of variables. In this setting, there are few relations in terms of essential variables, that 
is, variables that are in the support of the standard monomials associated to an ideal. The remaining ones are 
of the type Xi — g where the leading term Xi is not an essential variable and the support of g has only essential 
variables. Therefore computation of a Grobner basis can be restricted to a proper subring of R containing only 
essential variables. The algorithm identifies these variables and computes relations of the first type using the 
BMA. The relations of type Xi — g are computed using standard linear algebra techniques. We have implemented 
the algorithm, which we call EssBM, in Macaulay 2. 

The paper is organized as follows. First we describe the EssBM algorithm. In Section 3, we provide the 
theoretical support for the algorithm and include a complexity analysis. In Section 4, we compare its performance 
to the BMA, as implemented in Macaulay 2. We conclude our paper with a discussion of future directions. 



Let R = fc[iEi, . . . , x n ] where k is a field, and y be a fixed term order on R. Consider a variety V C k n of points with 
multiplicity one and |V| = m < oo. Here we are primarily interested in finite fields, where these conditions will 
automatically be satisfied for all varieties. The goal of the EssBM algorithm is to construct the reduced Grobner 
basis G with respect to y for the ideal I(V) of points in V and the set 23(G) of standard monomials associated 
to G, which forms a basis for the fc-vector space R/I(V). The algorithm constructs a set EV C {x\ . . . ,x n } of 
essential variables, a set SM of monomials on {x± . . . ,x n }, and subsets GB and Rel of the ring R. We will see 
below that G will be given by GB U Rel and 23(G) by the set SM. The support (defined in the next section) 
of the elements in SM is the set EV . We let EVi, SMi, GBi, and Reh denote the i-th approximations of the 
corresponding sets. 

Initialize each set as follows: EVo = {}, SMo = {la}, GBq = {}, and Relo — {}. Let [n] denote the set 
{1, . . . , n} and x a the monomial xi 1 ■ ■ ■ x^ n . For each i £ [n], do the following. Find the i-th smallest variable, 
say Xi. Suppose there are r monomials x ai , . . . , x ar in SMi-i that are smaller than Xi in the given ordering. Try 
to write Xi as a fc-linear combination of these monomials. That is, find (if they exist) ci, . . . , c r € k, where 



2 The EssBM Algorithm 
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and x a (t) is the evaluation of x a at the £-th point in V for t 6 [m]. If there are such coefficients, then 

r 

a*(«)-X>ir 0l (t) = u 

3=1 

for every t £ [m] and it follows that h := x t ~ Ylj=i c j xaj 6 I(^) n fc[-EVj-i U {xj}], where fc[_ET4_i U {2^}] is 
the polynomial ring in the variables in EVi-i U {a:;}. Since the monomials x aj were chosen so that Xi y x a ' ' , 
it follows that Xi is the leading term of an element of I(V) and so is not a standard monomial. In this case let 
Rek — Rek-i U {h}. If there is no solution to the system in {T}, then Xi is a standard monomial. In this case let 
EVi = EVi-i U {xi}, and compute the Grobner basis GBi and the set SMi of standard monomials for the ideal 
I(V) n k[EVi] of the points projected onto the variables in EVi. When i — n, return the sets G := GB n U Rel n 
and 58(G) := SM n . 

Below we give pseudo-code for the complete algorithm, which has been implemented in Macaulay 2. While the 
BMA computes separators for the points in V in addition to the Grobner basis and the set of standard monomials, 
the implementation in Macaulay 2 does not. In order to appropriately compare the two implementations, we do 
not include separators in this version of EssBM. However, our algorithm can easily be modified to return the 
separators at an additional cost of 0(m). 

For simplicity, let \xj{t)]YL\ denote the (m x l)-column vector 

/ ^(1) \ 



The EssBM Algorithm 

Input: V a variety; >^ a term order 

Output: G the reduced Grobner basis for I(V) with respect to 
23(G) the set of standard monomials for G 



1. Initialize: EVo := {}; SM := {1 R }; GB := {}; Rel := {}. 

2. For i from 1 to n do 

3. Xi := i-th smallest variable 

4. S := k[EVi-\ U {xi}] with term order induced by >^ 

5. r := |SMj_i| and LMi := {x aj Xi : x aj € SMi-i, 1 < j < r} the standard monomials less than Xi 

6. At := (m x (s + l))-matrix with s = \EVi-i\ 

first column [xi{t)]YL\ and s columns [xj (t)]t2=i for all Xj G EVi-i 

7. Evah := (m x r)-matrix (x aj (pt)), where x aj £ LMi is evaluated on p t , the point in row t of Ai 

8. If there is a solution c = (ci, . . . , c r ) T to the system of linear equations Evah • c = [^(t)]^! 

9. then i?e?i := i?e/i_i U {a;^ - J] Cja; a 3 } where x a ' G LMi 

10. else EVi ■= EVi-i U {x^} and compute GBi and SAL in k[EVi\ using the BMA on Ai 

11. Return G = GB„ U J?e«„ and 03(G) = SAf„ 

The variables in EV n are called essential. The polynomial Kj — X] c 3 :r£I; ' computed in the i-th step of the 
algorithm has Xi as its leading term since the monomials x aj were chosen to be smaller than Xi. The variables x% 
are called inessential since they can be written in terms of essential variables. 
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3 Theoretical Background 



In this section, we provide a detailed proof of the correctness and worst-case time complexity of the EssBM 
algorithm. Before stating and proving the main results, namely Theorems O [7] and [5] we begin with some 
preliminaries. 

Recall that the matrix Ai has rows corresponding to the points in V projected onto the coordinates defined 
by EVi = EVi-i U {xi}. Let Pi be this set of projected points. 

For the remainder of this paper, we use the shorthand notation I for the ideal I(V) and k[EVi] for the 
polynomial ring in the variables in the set EVi. Also, we let G — GB n U Rel n and 25(G) the set of standard 
monomials for G. 

Lemma 1. The equality l{Pi) = I D k[EVi\ holds. 

Proof. This follows immediately from the construction of the ideal T-(Pi). □ 

Corollary 2. The set GBi is the reduced Grobner basis for the ideal I n k[EVi] with respect to y and SMi is the 
set of standard monomials for I n fcf-EVi] with respect to GBi. In particular, the statement holds for i — n. 

Proof. The sets GBi and SMi are the reduced Grobner basis and the set of standard monomials, respectively, for 
the ideal I(Pi) in k[EVi\. From the previous lemma, we have that I(Pi) = Ink[EVi\. Hence the result follows. □ 

Let / £ R be a polynomial. We define the support of /, denoted by supp(f), to be the set of variables that 
appear in /. By construction, supp(f) is the smallest set X C {x\, . . . ,x n } such that / £ k[X]. The support of a 
set of polynomials S is the union over the support of each polynomial g £ S. Let LT(f) denote the leading term 
of / with respect to a given term order. The tail of / is the polynomial tail(f) := / — LT(f). 

Lemma 3. Let f £ R be such that supp(f) C EV n U {xp 1} . . . ,xp a } where xp 1 -<■■■-< Xp 3 are inessential 
variables. Suppose that supp(LT(f)) C EV„. Then there is f* £ R such that supp(f*) C EV„U {xp 1 , ... ,xp s _ 1 }, 
the polynomial /* has the same leading term as f, and f — f* 6 I. 

Proof. Consider the largest inessential variable xp, . We can write 

r 

/ = LT(/) + 

where supp(hi) C EV n U {x/3 lt . . . ,x/s s l }. As xp B is an inessential variable, there is an element xp s + g of Rel n 
with leading term xg a . Note that supp(g) C EV„. Define the polynomial /' from / by replacing each (x/3 a Y 
with -(xf 3s y- 1 g: 

r 

f = LT(f)-J2M i ~ 1 9h i . 

i=0 

Then 

r 
i=0 

since (Xf,,y + (x f) ,) i - 1 g=(xf),y- 1 (xi3,+g)el. As LT{f) y xp s y LT(g), we have that LT(f') = LT(f). Let/* 
be the polynomial obtained after r replacements of xp B . Note that we have / — /* 6 I and LT(f*) — LT(f). 
Since we have replaced all occurrences of xp e , it follows that supp(f*) C EV n U {xg 1 , ■ ■ ■ , xp B _ 1 }. □ 

This lemma gives us a way of removing inessential variables from a polynomial in / without affecting its leading 
term, which will be useful for proving the correctness of EssBM (Theorems [5] and [7)l . In fact, we can remove all 
inessential variables. We emphasize this fact with the following corollary. 

Corollary 4. Let f € R. Then there is f* £ R such that supp(f) C EV n , LT(f) = LT(f), and f - f e I. 
Theorem 5. The set G is the reduced Grobner basis for I with respect to y. 

Proof. We first show that G C I. Consider g £ G. If g £ GB n , then g £ I. Suppose that g £ Rel n . Then g is 
of the form Xi — CjX aj for some Cj £ k and x aj £ R — k[xi, . . . , x n ]. The coefficients Cj were chosen so that 
Xiit) = ^2cjX aj {t) for all t £ [m]. Therefore by construction g £ I. 

Now let / £ I. We must show that there is some g £ G such that LT(g) | LT(f). We distinguish two cases. 

Case 1: supp(LT{f)) (£ EV n . 

Suppose that LT(f) contains an inessential variable Xi. By construction of the set Rel n , there is an element g 
of Rel n C G with leading term Xi. It follows that LT(g) divides the leading term of /. 
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Case 2: supp{LT{f)) C EV n . 

Recall that the set GB n is a Grobner basis of the projection of I onto the variables in EV„ (see Corollary [5| ■ 
If supp(tail(f)) is also contained in EVn, then / G k[EV n ] and there is a g G GB n C G whose leading term 
divides LT(f). 

Assume that supp(tail(f)) <f_ EV„. Using Corollary [3] we can find ft G / such that supp(f — ft) C EVn 
and LT(f — ft) — LT(f). Since / — h €. k[EV n ], there is a g G GB„ C G whose leading term divides 
LT{f-h)=LT{f). 

To prove that G is reduced, let p ^ ft G G. We wish to show that <? and ft satisfy the following criterion: 

LT(g) does not divide any monomial in ft. (2) 
We consider the following four cases. 
Case 1: g,h G GB n . 

As G-B„ is the reduced Grobner basis for the ideal I projected onto the essential variables, then g, ft satisfy (2). 
Case 2: <?, ft G Rel n . 

Let LT(g) = Xj and h — Xj ~ £\ CiX a; for i 7^ j. Note that supp(h) C EVj—i U {xj}. Clearly a;* does not 
divide Xj. As supp(tail(h)) C .EVn and aii ^ BVn, then Xi does not divide any monomials in tail(h). 
Case 3: g G GB„ and ft G 7?ei n . 

Let LT(ft) = Xi for some inessential variable. This will not be divisible by LT(g), which contains at least one 
essential variable. All other terms x a of ft are standard monomials for the projection of / onto the variables 
in EVi\ in particular, supp(x a ) C EVi. It follows that if supp(g) C EVi, then LT(g) does not divide any 
term of ft. By Corollary [2] supp(g) contains only essential variables. Thus if supp(g) is not contained in EVi, 
then supp(g) must contain a variable Xj with Xi -< ay. This divides some term x of <?, and it follows that 
if LT(g) divides some term x a of ft, then Xj ^ x b < LM(g) X x a X Kj, which contradicts the assumption 
that Xi -< Xj . 

Case 4: j £ -Rei n and ft G GB n . 

Then LT(g) is some inessential variable, say Xi. However, supp(h) C -EV n and so p,ft satisfy criterion p). 

□ 

Next we compute the number of elements in 25(G) and show the relationship between 25(G) and the set SM n - 

Lemma 6. The set 23(G) has \V\ elements. 

The previous lemma is usually stated for algebraically closed fields k and proved with the help of the Strong 
Hilbert Nullstellensatz (see [H]). We include a proof of the statement for the case where all points have multiplicity 
one, as is being assumed throughout the paper. 

Proof. Suppose V = {ai, . . . ,a m } and define U := I({di}). Then / = I(U^Li{ a i}) = Dlii since each point m 
has multiplicity one. Note that each of the ideals Ii is maximal and it follows that they are pairwise comaximal. 
Consider the quotient ring R/I. By the Chinese Remainder Theorem, there is a ring homomorphism such that 

R/I R/h x ■ • ■ x R/Im- 

As each Ii is maximal, then each R/h = k and it follows that R/I = k m , as rings. Further, the quotient ring 
and fc m can be viewed as fe-vector spaces, and the isomorphism can be extended to an isomorphism of vector 
spaces. Hence, the dimension of R/I as a vector space is dirrik(R/ 1) = m. Since 25(G) forms a basis for the vector 
space R/I (Proposition 2.1.6 in [2]), we conclude that |25(G)| = m — |V|. □ 

Theorem 7. The set SM n is the set of standard monomials for I with respect to G. 

Proof. By Corollary [2] we have that SM n is the set of standard monomials for the ideal Ink[EV n ] with respect to 
the Grobner basis GB n . As V has finitely many points, then j23(G)| = \V\. Consider a monomial x a G 23(G). If 
x a £ k[EV n ], then it contains an inessential variable, say x^. As x, ; is the leading term of an element in Rel n C /, 
it is not a standard monomial for G, contradicting the assumption that x a 23(G). Therefore x a G fc[SVn]. 
By construction, x a LT(I). Using the set-containment relation 

LT(Ink[EV n ]) C LT(I), 

it follows that x a £ LT(Ink[EV n ]) and so 23(G) C SM„. To see equality, note that the set P n of projected points 
defined by EV n has at most as many points as V. Then ISAfnl = \Pn\ < |V| = m. Since 23(G) C SM„, it follows 
that m = |23(G)| < \SM n \ < m. Hence 23(G) = SM n ; that is, SM n is the set of standard monomials for I with 
respect to G. □ 
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We conclude this section with a complexity analysis of EssBM. 

Theorem 8. The EssBM algorithm terminates and has worst-time complexity 0(nm 3 ) + 0(m 6 ), which is domi- 
nated by 0(nm 3 ) when m <C n. 

Proof. We compute the complexity of each step and then provide a summary at the end. Step 1 has complexity 
O(l). In Step 2, the algorithm enters a loop of length n. Steps 3-8 are executed in each iteration of the loop. 
They have the following complexities: 

Step 3. O(l): Executing this step requires constant time since the variable order, given as part of the declaration of 
the term order, is maintained in one array. 

Step 4. 0(m 2 ): This step may not even be required by all implementations; if required, it involves passing 0(m 2 ) 
variables to a new object of size 0(m 2 ). 

Step 5. 0(m 3 ): As term orders are typically stored as matrices, in this case the term order -<s is a matrix of 
dimension 0(m 2 ). Determining the order between two monomials of S requires multiplication of a vector of 
length O(m) by this matrix. So for each monomial x a 6 SMi-i, there are at most m 2 operations required 
for comparing x% to x a and there are at most m such monomials. 

Step 6. 0(m): An (m x l)-column vector is added to a matrix with columns corresponding to the variables in EVi—i. 

Step 7. 0(m 3 ): As there are at most m variables in each monomial and at most m 2 entries in the matrix, the cost 
of executing this step is 0(m 3 ). 

Step 8. 0(m 3 ): Solving a linear system of m equations in r < m unknowns requires 0(m 3 ) time. 

Step 9 has complexity O(l) and will be executed at most n times. 

Since there can be at most m essential variables, Step 10 will be executed at most m times. The complexity 
of each execution of Step 10 is 0(m 5 ): Updating EVi is a constant operation. However, computing GBi and SMi 
for the matrix A% is associated to the cost of calling the BMA, which is quadratic in the number of variables and 
cubic in the number of points. In this case, the numbers of variables and points are given by the dimensions of Ai. 
Since both row and column dimensions are bounded above by m, it follows that the complexity of executing this 
step is 0(m 5 ). 

Step 11 has complexity 0(n + m 2 ): Note that there are 0(m 2 ) elements in GB n (see [50]), 0(n) relations in 
Rel n , and m monomials in SM n . So returning these sets requires 0(n + m 2 + m) operations. 
Hence, we can calculate the total complexity C(EssBM) of the algorithm as follows: 

C(EssBM) = O(l) + 0(n) [0(1 + rn 2 + rn + rn + m 3 + m 3 + 1)] + 0(m)0(m 5 ) + 0(n + m 2 + m) 
= 0(nm 3 ) + 0(m 6 ). 

When m <C n, then 0(nrn 3 ) becomes the dominating term and the above estimate reduces to 

C(EssBM) = 0(nm 3 ). 

□ 

4 Performance of the EssBM Algorithm 

To test the performance of our algorithm, we compared its run-time to that of the as implemented in 

Macaulay 2, on randomly generated varieties in k n . For this analysis, we let the field k be F p for p £ {3, 17}. 
Since the complexities of the two algorithms depend on m and n, we chose a range of values for these parameters, 
namely, m £ {5, 10, 15} and n £ {100, 150, 200, 250, 300}. For each set of parameters p, m, and n, we generated 
10 varieties using a built-in random number generator in Macaulay 2, without specifying prior constraints on the 
relative position of the points in the variety. We performed this experiment using two term orders: a lexicographic 
order (lex) and a graded reverse lexicographic order (grevlex), each with the same variable order. 

Figures [T] and [5] show the run-times for the two algorithms for p = 3 and m = 5, 15. As the run-times for 
m — 10 fall between the m = 5 and m = 15 settings, we omitted them from the plots. We display the results for 
all parameters settings in the appendix. The run-times for p = 17 are similar. 

As a measure of the stability of the run-time data, we computed the coefficient of variation, defined to be the 
ratio of the standard deviation to the mean of the data. For the grevlex experiments, this coefficient ranges from 

lr The Buchberger-Mollcr algorithm has been implemented as the function points in the "Points" package of Macaulay 2 distribution 
version 0.9.8. 
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0.004 to 0.2, whereas for the lex experiments it ranges from 0.01 to 0.1. Since this implies very low variability of 
the run-times for fixed p, n, and m, we displayed only mean values in Figures Q] and [5] 

The empirical results corroborate our theoretical prediction that for m Cn, the EssBM algorithm outperforms 
the BMA. For small n, however, we observe that EssBM is slower, which we attribute to the overhead costs 
associated to multiple calls to the BMA. 

5 Discussion 

Recently, applications of Grobner bases as a promising model selection tool in molecular biology have been 
proposed 3, 18]. These applications require computation of a Grobner basis for a zero-dimensional ideal I(V) 
in a polynomial ring k[xi,... ,x n ], where \V\ — m <C n. Previously, no algorithms for computing Grobner bases 
optimized for m <^ n had been available. The run-time of the existing implementations was a bottleneck in 
applications of the methods in [3] and |18| to data sets whose size is of the order typical for biochemical data sets 
such as microarray data. 

The EssBM algorithm presented here goes some way towards alleviating this problem in that it reduces the 
worst-case complexity, which is 0(n 2 m 3 ) for the standard Buchberger-Moller algorithm, to 0(nm 3 ) for m <C n. 
Our implementation and testing indicate that for a small number of distinct points in general position, EssBM 
starts outperforming a standard implementation of the BMA when the number of variables exceeds 200. This 
should make it possible to use the methods of [3] and [18] for analysis of larger data sets than was hitherto possible. 
Unfortunately, the worst-time complexity estimate 0{nrn i + m 6 ) of the EssBM algorithm suggests that it may 
still be infeasible for moderately large m. We are currently working on a related algorithm that would further 
reduce this complexity. 
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Figure 4: Run-times for 10 randomly generated varieties and grevlex. 
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